A_GTP_exploring_phenotype_data

.libPaths( c( "/opt/R/3.6.2/lib/R/library" , .libPaths() ) )
.libPaths()

[1] "/opt/R/3.6.2/lib/R/library"                      
[2] "/PHShome/bm363/R/x86_64-pc-linux-gnu-library/3.6"

load datasets

library(tidyverse)
library(readxl)

GTP_pheno3 <- read_excel("/data/humgen/GTP_genomics/PLINK/input/pheno/ResslerData_cleaned_reduced_For Machine Learning_withSIDs.xlsx", sheet = 1) %>% 
  print() 

GTP_pheno3 %>% 
#  count(sex)
#  count(PTSD_diagnosis)
  subset(select=c("sex","PTSD_diagnosis")) %>% table()

GTP_pheno3 <- read_excel("/data/humgen/GTP_genomics/PLINK/input/pheno/ResslerData_cleaned_reduced_For Machine Learning_withSIDs.xlsx", sheet = 1) %>% 
  print() 

GTP_pheno3 %>% 
#  count(sex)
#  count(PTSD_diagnosis)
  subset(select=c("sex","PTSD_diagnosis")) %>% table()

      PTSD_diagnosis
sex       0    1 2092 7011
  0    1300  552    0    0
  1    3517 1631    0    0
  519     0    0    1    0
  8584    0    0    0    1

What is a better trauma exposure measure we should use? Do we have CRP and other molecular markers measured? Do we have BMI and other anthropometric measurements?

GTP_pheno4 <- GTP_pheno3 %>% 
#  filter(!is.na(PTSD_status)) %>% 
#  filter(PTSD_status<125) %>% 
  subset(select=c("PTSD_diagnosis", "CTQTOT","education")) %>% 
  mutate(CA4 = cut(CTQTOT, breaks = c(-Inf,40,60,80,Inf), right=T)) %>%
  mutate(EA4 = cut(education, breaks = c(-Inf,1,3,4,Inf), right=T)) 

GTP_pheno4 %>% print()
#, labels=c("low","midlow","midhigh","high")

GTP_pheno4 <- GTP_pheno3 %>% 
#  filter(!is.na(PTSD_status)) %>% 
#  filter(PTSD_status<125) %>% 
  subset(select=c("PTSD_diagnosis", "CTQTOT","education")) %>% 
  mutate(CA4 = cut(CTQTOT, breaks = c(-Inf,40,60,80,Inf), right=T)) %>%
  mutate(EA4 = cut(education, breaks = c(-Inf,1,3,4,Inf), right=T)) 

GTP_pheno4 %>% print()

#, labels=c("low","midlow","midhigh","high")

Genotype data

# Number of variants and samples
wc -l /data/humgen/GTP_genomics/PLINK/input/PGC_GTP_June2017/GTPC/qc/pts_gtpc_mix_am-qc.{bim,fam}

# an alternative geno data file location (befoe QC with) with slightly more samples and variants
wc -l /data/humgen/GTP_genomics/PLINK/input/PGC_GTP_June2017/GTPC/starting_data/Adam_final_GTP_pheno.{bim,fam}

GTP_fam1 <- read_delim("/data/humgen/GTP_genomics/PLINK/input/PGC_GTP_June2017/GTPC/qc/pts_gtpc_mix_am-qc.fam",
                         col_names = F, delim = "\t") %>% 
  rename_all(funs(c("FID","IID","p1","p2","gend","phen"))) %>% 
  mutate(ID = str_replace(IID,"PTSD_",""))

GTP_pheno4 <- GTP_pheno3 %>% 
  mutate(sid = as.character(sid)) %>% 
  inner_join(GTP_fam1, by=c("sid"="ID")) %>% 
  print()

#write.table(GTP_pheno4, "/data/humgen/burook/GTP_analysis/GTP_pheno1.txt", row.names = FALSE, col.names = TRUE, quote = FALSE, sep = "\t")

Ancestry estimation and confirmation

Prepare a file with self-identified race.

GTP_fam1 <- read_delim("/data/humgen/GTP_genomics/PLINK/input/PGC_GTP_June2017/GTPC/qc/pts_gtpc_mix_am-qc.fam",
                         col_names = F, delim = "\t") %>% 
  rename_all(funs(c("FID","IID","p1","p2","gend","phen"))) %>% 
  mutate(ID = str_replace(IID,"PTSD_",""))

Parsed with column specification:
cols(
  X1 = col_character(),
  X2 = col_character(),
  X3 = col_double(),
  X4 = col_double(),
  X5 = col_double(),
  X6 = col_double()
)

GTP_pheno4 <- GTP_pheno3 %>% 
  mutate(sid = as.character(sid)) %>% 
  inner_join(GTP_fam1, by=c("sid"="ID")) %>% 
  print()

#write.table(GTP_pheno4, "/data/humgen/burook/GTP_analysis/GTP_pheno1.txt", row.names = FALSE, col.names = TRUE, quote = FALSE, sep = "\t")

GTP_race2 <- GTP_pheno4 %>% 
  subset(select=c("IID","race_ethnic")) %>% 
  rename(Ethnicity_race=race_ethnic)

#write.table(GTP_race2, "/data/humgen/burook/GTP_analysis/GTP_race2.txt", row.names = FALSE, col.names = TRUE, quote = FALSE, sep = "\t")

t1 <- read_delim("/data/humgen/burook/GTP_analysis/table_pred", col_names = TRUE, delim = " ")

t1 %>% 
  filter(is.na(SuperPopulation)) %>% 
  group_by(predtree) %>%
  summarise(n=n()) %>%
  mutate(n_perc = 100*n / sum(n)) %>% 
  print()

Computing PRS

 t1 %>% 
   left_join(target_race, by=c("Sample"="IID")) %>% 
   filter(!is.na(Ethnicity_race.y)) %>% 
   subset(select=c("Ethnicity_race.y","predtree")) %>% 
   table() %>% 
   print()

                predtree
Ethnicity_race.y  AFR  EAS  EUR  SAS
            0    5698    0   30    4
            1      13    0   17   12
            2       0    1    3    4
            3       1    0  191    2
            4      82    0   19    9
            5      67    0    9    0
            8504    2    0    0    0

# with freeze-2 PGC-PTSD summ stats
R --file=/PHShome/bm363/bin/PRSice_v1.25/PRSice_v1.25.R -q --args  \
  plink /PHShome/bm363/bin/plink2 \
  base /data/humgen/burook/sum_stats/pts_aam_freeze2_overall.results \
  target /data/humgen/burook/GTP_analysis/gtp_ptsd/tmp_gtp1 \
  pheno.file /data/humgen/burook/GTP_analysis/GTP_race2.txt \
  binary.target T \
  covariates C1,C2,C3,C4,C5,C6,C7,C8,C9,C10 \
  clump.snps T \
  prune.snps F \
  quantiles T \
  figname gtp_ptsd \
  wd /data/humgen/burook/GTP_analysis/gtp_ptsd/

rm /data/humgen/burook/GTP_analysis/gtp_ptsd/rawfile.raw
ls /data/humgen/burook/GTP_analysis/gtp_ptsd/
#head /data/humgen/burook/GTP_analysis/GTP_pheno1.txt